#ifndef LINEAGE_H
#define LINEAGE_H

#include <math.h>
#include <vector>
#include <string>
#include <iostream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
#include "tree.h"
#include <time.h>

#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf.h>
#include <gsl/gsl_sf_gamma.h>

gsl_rng * r2;
//long seed;
//seed = rand()*time (NULL) * getpid();    
//gsl_rng_set(r, seed);   

// lineage class contains matrices for data, lineage structure, data, pool structure, cones and 
// vectors for cones, cone locations and the random nodes used in proposing a tree.  it also 
// includes a tree, the number of lineages, pools and SNPs.
// functions for reading data, setting variables, initializing the structures and normalizing vectors

class lineage
{
	public:
//		basic objects: tree, pools, lineages, cones
		tree treeT;
		vector < vector <double> > poolS;	//pool structure
		vector < vector <double> > lineageR;	//lineage structures

		int numLineage;				// number of lineages
		int numPools;				// number of pools
		int numSNPs;				// number of SNPs

		double logLineageProposal;
		double xi;
		double lambda; // the Poisson parameter
		double eta; // probability that a SNP is poisson  
	
// 		constructors/deconstructors/operators
		lineage(){lineageR.clear(); poolS.clear(); 
			  numLineage = 0; numPools = 0; numSNPs = 0; 
			  r2 = gsl_rng_alloc(gsl_rng_rand48); long seed; seed = rand()*time (NULL) * getpid();  gsl_rng_set (r2, seed); }
		~lineage(){};
		lineage& operator=(const lineage&);
	
		//	lineage(vector< vector <double> > S,vector< vector <double> > D);
		//	DataIO;
		void readData(string data_file, vector< vector< vector <int> > > & dataD);	//reads data
		void setNumLineages(int k);		//sets numLineage
		void setNumPools(int);			//sets numPools
		void setNumSNPs(int);			//sets numSNPs
		void setXi(double);
		void setEta(double);
		void setLambda(double);

		// Initialization Routines
		void initializeLineages(void);	//generates random lineageR for initialization
		void initializePools(void);		//initializes the pools randomly	
		void initializeParameters(const vector< vector< vector <int> > > & dataD);	//initializes all the parameters, once data and numLineage 

		// Model refinements and math functions
		void normalizeModel(void);		//performs normalization so tree sum to 1
		double logBinomial(long i, long j);	//finds the log of binomial
		double logFactorial(long);		//finds the log factorial, using the stiriling approx for large values
		void normalizeVector(vector<double> &); //normalizes a vector

		//Proposal routines
		void proposePool(int i);		// also proposes a pool, but stupidly
		void proposeLineage(double,vector< vector< vector <int> > >&);
		void proposeLineage(vector< vector< vector <int> > >& dataD);
		void proposeEta();
		void proposeLambda();
		void proposeTree();
		void proposeXi();
		double calculateSiteLogLikelihood(vector< vector< vector <int> > >&,int,vector<double>);
		double calculateSiteLogLikelihood(vector< vector< vector <int> > >&,int);
		double calculateSiteLogBinomial(int, int);

		void proposePool(vector< vector< vector <int> > >& dataD, int i, int iter);
		vector<double> proposeSiteLineage(vector< vector< vector <int> > >&,int,int);

		// Useful subroutines + output
//		void createLineageGraphHomoplasy(string filename);
		int vectorRandomElement(vector<double>);
		void calculateLogProposal(vector< vector< vector <int> > >& dataD);
		string createLineageGraph();
		double modelLengthBelowNode(int j);
		
	private:
		vector<double> generateRandomState(int);
		double glog(double);
		vector < double > convertBase(int,int,int);
		double matrixExp(double t, int start, int end);
		double recursionProb(int node, int j, int state);
};


//overloaded output
//output of lineage class to check for consistency
ostream& operator << (ostream& os, const lineage& s)
{  
	cout<<"\t Tree has "<<s.treeT.treeNodes.size()<<" nodes."<<endl
	    <<"\t Lineage: "<<s.lineageR.size()<<" "<<s.lineageR[0].size()<<endl
	    <<"\t Pool: "<<s.poolS.size()<<" "<<s.poolS[0].size()<<endl;
	int i,j,k,N = s.numPools, K = s.numLineage, M = s.numSNPs;
	
	double total;
	cout<<endl;
	cout<<"P:"<<endl;
	for(i = 0;i<N;i++)
	{
		total = 0;
		for(k=0;k<K;k++)
		{
			cout<<s.poolS[i][k]<<" ";
			total = total+s.poolS[i][k];
		}
		cout<<total;
		cout<<endl;
	}
	cout<<"R"<<endl;
	for(j =0; j<10; j++)
	{
		for (k = 0; k<K;k++)
		{	cout<<s.lineageR[k][j]<<" ";	}
		cout<<endl;
	}
	cout<<"Xi:"<<s.xi<<endl;
}

//overloaded assignment operator
lineage& lineage::operator=(const lineage &rhs) 
{
    // Check for self-assignment!
    if (this == &rhs)      // Same object?		void updateLineage();
      return *this;        // Yes, so skip assignment, and just return *this.
    else
    {
	treeT    = rhs.treeT;
    	poolS    = rhs.poolS;
	lineageR = rhs.lineageR;

	numLineage = rhs.numLineage;
	numPools   = rhs.numPools;
	numSNPs    = rhs.numSNPs;

	xi	 = rhs.xi;
	lambda   = rhs.lambda;
	eta      = rhs.eta;

	logLineageProposal = rhs.logLineageProposal;	

    	return *this;
     }
}

// normalizes tree and cone length to sum to 1
void lineage::normalizeModel(void)
{
	treeT.makeTreeUltrameric();
}

void lineage::setXi(double value)
{
	xi = value;
}

//sets the number of SNPs
void lineage::setNumSNPs(int num)
{
	numSNPs = num;
}

void lineage::setEta(double value)
{
	eta = value;
}

void lineage::setLambda(double value)
{
	lambda = value;
}

// fucntion initializes the lineage, pool, cones and tree with random values
// checks to see if data is entered and numLineage >0

void lineage::initializeParameters(const vector< vector< vector <int> > > & dataD)
{
//	cout<<"Data has size: "<<dataD.size()<<endl;
	if ((dataD.size()>0)&&(numLineage>0))
	{
		setNumSNPs(dataD.size());
//		setNumSNPs(10);
		setNumPools(dataD[0].size());
	
		cout<<"Sizes set..."<<endl;
		treeT.generateRandomTree(numLineage);
		treeT.hashLN();

		cout<<"Tree initialized.."<<endl;

	//	normalizeModel();	
		cout<<"Random tree generated.  Here it is:"<<endl;		
		treeT.displayTree();
		cout<<"Randomizing positions..."<<endl;
		initializePools();
		cout<<"Pools initialized..."<<endl;	
		initializeLineages();
		cout<<endl<<"Lineages initialized..."<<endl<<"Intialized."<<endl;
		setXi(1);
		setEta(0.01);
		setLambda(0.05);

 	}
	else
	{
		cerr<<" ** Some or all data not entered. \n ** All models have data. \n ** Otherwise they are just castles in the air."<<endl;
		exit(-1);
	}
}

//randomly initializes the pools, drawing from a dirichlet with a vector of 1s as parameter
void lineage::initializePools(void)
{
//	cout<<"Starting pools	test = gsl_ran_flat(r2,0,1);..."<<endl;
	int i,k;

	vector <double> dummy;
	dummy.clear();
	dummy.resize(numLineage,0);

	poolS.clear();
	poolS.resize(numPools,dummy);
	
	double alphaV[numLineage], theta[numLineage];
 
	for (i = 0; i<numLineage; i++)
	{
		alphaV[i] = 1;
	}

	for (i = 0; i < numPools; i++)
	{	
		gsl_ran_dirichlet(r2, numLineage, alphaV,theta);
		for (k = 0; k<numLineage;k++)
		{	
			dummy[k] = theta[k];
		}	
		normalizeVector(dummy);	
		for (k = 0; k < numLineage; k++)
		{
			poolS[i][k] = 1/(double)numLineage;//dummy[k];
		}
	}

}

void lineage::normalizeVector(vector<double> & vec)
{
	double total = 0;
	int i =0, check;

	for (i =0; i<vec.size(); i++)
	{	total += vec[i];	}
	
	if (total > 0)
	{
		for (i =0; i<(int)vec.size();i++)
		{
			vec[i] = vec[i]/total;
		}
	}
	else
	{
		for (i =0; i<(int)vec.size();i++)
		{
			vec[i] = 1/(double)vec.size();
		}		
	}
	check = 0;
	for (i = 0; i<vec.size();i++)
	{
		if (vec[i]==0)
		{	check = 1; vec[i] = 1e-20;	}
	}
	if (check == 1)
	{	normalizeVector(vec);	}
	//return vec;
}

void lineage::proposeEta()
{
	double test = gsl_ran_flat(r2,0,1), proposed_eta, current_eta = eta,value, max_val = 0.25,min_val=0.0;
	logLineageProposal = 0;

	logLineageProposal = glog(gsl_ran_beta_pdf(eta,1,100));
	eta = gsl_ran_beta(r2,1,100);
	logLineageProposal -= glog(gsl_ran_beta_pdf(eta,1,100));
/*

	if (test < 0.5)
	{
		// then propose below
		proposed_eta = current_eta - (current_eta-min_val)*gsl_ran_beta(r2,1,20);
		value = (current_eta-proposed_eta)/(max_val-proposed_eta);
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (current_eta-proposed_eta)/(current_eta-min_val);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
	else
	{
		proposed_eta = current_eta + (max_val - current_eta)*gsl_ran_beta(r2,1,20);
		value = (proposed_eta-current_eta)/(proposed_eta-min_val);
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (proposed_eta-current_eta)/(max_val-current_eta);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
*/
//	eta = proposed_eta;
}

void lineage::proposeXi()
{
	logLineageProposal = glog(gsl_ran_exponential_pdf(xi,1));;
	xi = gsl_ran_exponential(r2,1);
	logLineageProposal -= glog(gsl_ran_exponential_pdf(xi,1));;	

	double test = gsl_ran_flat(r2,0,1), proposed_xi, current_xi = xi, value;
	logLineageProposal = 0;


	
	if (test < 0.5)
	{
		// then propose below
		proposed_xi = current_xi-current_xi*gsl_ran_beta(r2,1,20);
		value = (current_xi-proposed_xi)/(proposed_xi);
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (current_xi-proposed_xi)/(current_xi);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
	else
	{
		proposed_xi = current_xi+(current_xi)*gsl_ran_beta(r2,1,20);
		value = (proposed_xi-current_xi)/proposed_xi;
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (proposed_xi-current_xi)/(current_xi);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
	xi = proposed_xi;

}


void lineage::proposeLambda()
{
	double test = gsl_ran_flat(r2,0,1), proposed_lambda, current_lambda = lambda,value;
	logLineageProposal = 0;

	logLineageProposal += glog(gsl_ran_beta_pdf(lambda,1,20));
	lambda = gsl_ran_beta(r2,1,20);
	logLineageProposal -= glog(gsl_ran_beta_pdf(lambda,1,20));

/*
	if (test < 0.5)
	{
		// then propose below
		proposed_lambda = current_lambda-current_lambda*gsl_ran_beta(r2,1,20);
		value = (current_lambda-proposed_lambda)/(1-proposed_lambda);
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (current_lambda-proposed_lambda)/(current_lambda);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
	else
	{
		proposed_lambda = current_lambda+(1-current_lambda)*gsl_ran_beta(r2,1,20);
		value = (proposed_lambda-current_lambda)/proposed_lambda;
		logLineageProposal += glog(gsl_ran_beta_pdf(value,1,20));
		value = (proposed_lambda-current_lambda)/(1-current_lambda);
		logLineageProposal -= glog(gsl_ran_beta_pdf(value,1,20));
	}
*/
//	lambda = proposed_lambda;
}

vector<double> lineage::generateRandomState(int K)
{
	
	int i;
	double test;
	vector <double> outData;
	
	outData.clear(); outData.resize(K,0);

	for (i = 0; i < K; i++)	
	{
		test = gsl_ran_flat(r2,0,1);
		if (test < 0.5)
		{	outData[i] = 1.0;	}
	}
	return(outData);

}


double lineage::calculateSiteLogBinomial(int a , int b)
{
	int dummy = 0;
	if (b < a){ dummy = a; a = b; b = dummy;	}
	double value = glog(gsl_ran_binomial_pdf(a, lambda, a+b));
	return(value);
}



double lineage::recursionProb(int node, int j, int state)
{
	int k = 0, r_node, l_node, i;
	double r_branch, l_branch,value = 0, value_l=0, value_r=0;

	if (treeT.isTip(node))
	{
		k = treeT.lineageNodeHash[node];
		if (state == lineageR[k][j])
		{	value = 1;	}
		else	
		{	value = 0;	}
	}
	else
	{

		r_node = treeT.treeNodes[node].right;
		l_node = treeT.treeNodes[node].left;

		r_branch = treeT.treeNodes[node].rightBranch;
		l_branch = treeT.treeNodes[node].leftBranch;

		value_l = 0;
		value_r = 0;
		for (i = 0; i < 2; i++)
		{

			value_r += matrixExp(r_branch,state,i)*recursionProb(r_node,j,i);
			value_l += matrixExp(l_branch,state,i)*recursionProb(l_node,j,i);

		}
		value = value_r*value_l;
	}
	return(value);
}

double lineage::matrixExp(double t, int start, int end)
{
	double value = 0;
	if (start == end)
	{	value = 0.5 + 0.5*exp(-2*xi*t);	}
	else
	{	value = 0.5 - 0.5*exp(-2*xi*t);	}
	return(value);
}

vector < double > lineage::convertBase(int number, int base, int length)
{
	int k;
	vector < double > out;
	out.clear(); out.resize(length,0);

	for (k = length; k > 0; k--)
	{
		if (pow(base,k) > number)
		{
			out[k-1] = 0;
		}
		else
		{
			out[k-1] = 1;
			number = number - pow(base,k);
		}
	}
	return(out);
}

void lineage::proposeLineage(double some, vector< vector< vector <int> > >& dataD)
{
	int k, K=numLineage, j, M = numSNPs,check = 1;
	double test = 0;
	vector <double> dummy,store;

	j = floor(gsl_ran_flat(r2,0,M-epsilon));
	k = floor(gsl_ran_flat(r2,0,K-epsilon));

	logLineageProposal = 0; 

	lineageR[k][j] = fabs(1-lineageR[k][j]);
}

void lineage::proposePool(int i)
{
	vector<double> outputVector;

	outputVector.clear();
	outputVector.resize(numLineage,0);

	int k =0, cand1 = 0, cand2 = 0, check;
	double uni = 0, combinedMass = 0;

	cand1 = floor(gsl_ran_flat(r2,0,numLineage-epsilon));
	check = 1;

	while(check)
	{
		cand2 = floor(gsl_ran_flat(r,0,numLineage-epsilon));
		if (cand1 !=cand2)
		{	check = 0;	}	
	}

	combinedMass = poolS[i][cand1]+poolS[i][cand2];

	uni = gsl_ran_flat(r2,epsilon,1-epsilon);	
	for (k = 0; k<numLineage;k++)
	{
		if (cand1 ==k )
		{	outputVector[k]  = combinedMass*uni;}
		else
		{
			if (cand2 == k )
			{	outputVector[k]  = combinedMass*(1-uni);}
			else
			{	outputVector[k] = poolS[i][k];	}
		}
		
	}			
	normalizeVector(outputVector);

	for (k = 0;k<numLineage;k++)
	{
		poolS[i][k] = outputVector[k];
	}

}

void lineage::proposePool(vector< vector< vector <int> > >& dataD, int i,int iter)
{
	int M = numSNPs, j=0, K = numLineage, N = numPools,  k;
	double count = 0, total = 0,tune = 0.1;
	
	vector<double> dummy;
	dummy.clear();  
	
	double alphaV[K];
	double theta[K];

	logLineageProposal = 0;
	tune = 0.5*(iter < 20000)+0.1*(iter < 40000)+1*(iter<60000)+2*(iter > 100000);

	for(k =0;k<K;k++)
	{	
		alphaV[k] = 0;
		for (j = 0; j < M; j++)
		{		
		  alphaV[k] += (lineageR[k][j])*dataD[j][i][0] + (1-lineageR[k][j])*dataD[j][i][1];
		  total += (double)dataD[j][i][0] + (double)dataD[j][i][1];
		}

		alphaV[k] = (alphaV[k]/total)*tune+1;

	}

	gsl_ran_dirichlet(r2, K, alphaV, theta);

	for (k = 0; k<K; k++)
	{
		logLineageProposal += (alphaV[k]-1)*gsl_sf_log(poolS[i][k]/theta[k]);
	}

	for (k = 0; k<K; k++)
	{
		poolS[i][k] = theta[k];
	}

}

void lineage::proposeTree()
{
	treeT.randomSPR();
	treeT.hashLN();
}

int lineage::vectorRandomElement(vector<double> probVec)
{
	int i = 0, check = 1, counter = 0;
	double  uni;
	vector<double> sumVec;

	sumVec.clear();	
	sumVec.resize(probVec.size(),0);


	for (i = 0; i<probVec.size(); i++)
	{
		if (i > 0)
		{	sumVec[i] = probVec[i] + sumVec[i-1];}
		else
		{	sumVec[i] = probVec[i];	}
	}
	


	uni = gsl_ran_flat(r2,0,1 - 1e-5);

	while(check)
	{
		if (uni < sumVec[counter])
		{	check = 0;	}		
		else
		{	
			if (counter== probVec.size()-1)
			{	check = 0;	}
			else
			{	counter ++;	}
		}
	}

	return counter;
}


// initializes random lineage structure
// returns the randomNodes needed to randomize (may be unncessary)
void lineage::initializeLineages(void)
{
	int j, M = numSNPs, k, K =numLineage;
	double test;
	vector < double > dummy;


	dummy.clear();
	dummy.resize(M,0);

	lineageR.clear();
	lineageR.resize(K,dummy);

	for (j = 0; j < M; j++)
	{
		for(k = 0; k<K; k++)
		{
			test = gsl_ran_flat(r2,0,1);
//			cout<<test<<endl;
			if ( k > 0)
			{
				if ((j+1) %(k+1) ==0)
				{
					lineageR[k][j] = 1;
				}
			}
			else
			{
				lineageR[k][j] = k % 2;
			}
		}
	}

}

// sets the numLineages
void lineage::setNumLineages(int k)
{
	numLineage = k;
}

// reads the dataf
// very complicated and unreadable
// but works
// please don't mess with 
void lineage::readData(string data_file, vector< vector< vector <int> > >& dataD)
{
	ifstream dataStream;
	dataStream.open(data_file.c_str());

	char buffer[10000];
	string sample;
	int i, j,counter, check;

	dataD.clear();
	
	vector<int> dummyValues;
	vector<int> placeHolder;
	vector< vector<int> > dummyVector;

	
	dummyVector.clear();
	placeHolder.clear();
	placeHolder.resize(2,-1);

	cout<<"Reading data..."<<endl;
	if (dataStream.is_open())
	{
		while(dataStream.good())
		{	
			dummyValues.clear();
			dataStream.get(buffer,1e4,'\n');
			//cout<<buffer<<endl;
			if (dataStream.good() == 0)
			{	break;	}
			check = 1; counter = 0;i =0;
			while(check)
			{
				//cout<<buffer[counter]<<" "<<int(buffer[counter])<<endl;4
				if (((int(buffer[counter])<=57)&&(int(buffer[counter])>=48))||(int(buffer[counter]) ==44)||(int(buffer[counter]) ==13)||(int(buffer[counter]) ==0))
				{
					if (int(buffer[counter])==44)
					{

						sample = buffer[counter-i];
						for (j = 1; j<i;j++)
						{	
							sample = sample + buffer[counter-i+j];
				
						}
    				                dummyValues.push_back(atoi(sample.data()));
						i = 0;					
					}
					else
					{
						if ((int(buffer[counter])==13)||(int(buffer[counter]) ==0))
						{
							sample = buffer[counter-i];
							for (j = 1; j<i;j++)
							{	
								sample = sample + buffer[counter-i+j];
				
							}

							dummyValues.push_back(atoi(sample.data()));				
							check = 0;						
						}
						i = i + 1;
						
					}
				}
				else
				{		
					check = 0;
				}
				
				counter = counter + 1;

			 }

			dummyVector.resize(dummyValues.size()/2,placeHolder);			
			dataD.push_back(dummyVector);
			for (i =0; i<(int)dummyValues.size();i++)
			{
				if (((double)i) <((double)dummyValues.size()/2))
				{	
					dataD[dataD.size()-1][i][0] = dummyValues[i];			
				}
				else
				{
					dataD[dataD.size()-1][i-dummyValues.size()/2][1] = dummyValues[i];
				}
			}
			dataStream.get();
		}
		dataD.pop_back();
	}
	else
	{
		cout<<"No data file with that name"<<endl;
		exit(-1);
	}
	cout<<"Data read."<<endl;
	dataD.pop_back();
	dataStream.close();
	//return &dataD;
}

//sets the number of pools
void lineage::setNumPools(int num)
{
	numPools = num;
}


string lineage::createLineageGraph()
{
	stringstream ss;
	
//	vector <int> homoplasyCounts = getHomoplasyCounts();
		
	int i = 0, T = 2*numLineage-3, N = numPools, K = numLineage, k,s;

	ss<<"graph G {"<<endl;
	ss<<"{ graph [orientation = landscape, overlap = true] }"<<endl; 
	
	double factor_edge = 30, factor_node = 1, min_edge = 1.5, node_size =  0.33, total, pie_size = 1;
	
	treeT.hashLN();
	
	ss<<"{ node [ shapefile = \"legend.ps\"] o } "<<endl;

	for (i = 0; i<treeT.treeNodes.size(); i++)
	{
		if (treeT.isTip(i))
		{
			ss<<"{ node [fontsize = 0, height = 0.25, width = 0.25, shape = circle, color = white, image = \"stars_"<<treeT.lineageNodeHash[i]+1<<".ps\"] "
				    <<"s_"<<treeT.lineageNodeHash[i]+1<<" }"<<endl;
		}	
	}
	for( i = 0; i<treeT.treeNodes.size(); i++)
	{
		ss<<"{ node [shape = circle, style = filled";
		if (treeT.isTip(i))
		{
			ss<<", color = cadetblue4, width = 1, fontsize = 50] "<<treeT.lineageNodeHash[i]+1<<" }"<<endl;
		}
		else
		{
		    ss<<", color = dodgerblue4, width = 0.1] "<<"n_"<<i<<" }"<<endl;
		}
	}	
	for (i = 0; i <treeT.treeNodes.size(); i++)
	{
		if (i == 0)
		{	
			ss<<"{ edge [ label = "<<treeT.treeNodes[i].parent <<", fontsize =30,";	
			if (treeT.treeNodes[i].parentBranch*factor_edge > min_edge)
			{	ss<< " width = 4, color = darkslategray4, len = "<<treeT.treeNodes[i].parentBranch*factor_edge<<"] ";}
			else
			{	ss<<" width = 4, color = darkslategray4, len = "<<min_edge<<"] ";}

			if(treeT.isTip(treeT.treeNodes[i].parent))
			{      ss<<" n_"<<i<<" -- "<<treeT.lineageNodeHash[treeT.treeNodes[i].parent]+1<<"}"<<endl;}
			else
			{      ss<<" n_"<<i<<" -- n_"<<treeT.treeNodes[i].parent<<"}"<<endl;}
		}
		if(treeT.isTip(i))
		{	
			ss<<"{edge [len = 4, style = invis] "<<treeT.lineageNodeHash[i]+1<<" -- s_"<<treeT.lineageNodeHash[i]+1<<"}"<<endl;
		}
		else
		{
			if (treeT.treeNodes[i].rightBranch*factor_edge >min_edge)
			{
				ss<<"{ edge [ width = 4, color = royalblue4, len = "<<factor_edge*treeT.treeNodes[i].rightBranch<<", label = "<<treeT.treeNodes[i].right<<", fontsize =30]";	

				if(treeT.isTip(treeT.treeNodes[i].right))
				{      ss<<" n_"<<i<<" -- "<<treeT.lineageNodeHash[treeT.treeNodes[i].right]+1<<"}"<<endl;}
				else
				{      ss<<" n_"<<i<<" -- n_"<<treeT.treeNodes[i].right<<"}"<<endl;}					
			}
			else
			{
				ss<<"{ edge [ width = 4, color = royalblue4, len = "<<min_edge<<", label = "<<treeT.treeNodes[i].right<<", fontsize =30]";	

				if(treeT.isTip(treeT.treeNodes[i].right))
				{      ss<<" n_"<<i<<" -- "<<treeT.lineageNodeHash[treeT.treeNodes[i].right]+1<<"}"<<endl;}
				else
				{      ss<<" n_"<<i<<" -- n_"<<treeT.treeNodes[i].right<<"}"<<endl;}
			
			}
			if (treeT.treeNodes[i].leftBranch*factor_edge >min_edge)
			{	
				ss<<"{ edge [ width = 4, color = royalblue4, len = "<<factor_edge*treeT.treeNodes[i].leftBranch<<" , label = "<<treeT.treeNodes[i].left<<", fontsize =30]";	

				if(treeT.isTip(treeT.treeNodes[i].left))
				{      ss<<" n_"<<i<<" -- "<<treeT.lineageNodeHash[treeT.treeNodes[i].left]+1<<"}"<<endl;}
				else
				{      ss<<" n_"<<i<<" -- n_"<<treeT.treeNodes[i].left<<"}"<<endl;}
			}
			else
			{
				ss<<"{ edge [ width = 4, color = royalblue4, len = "<<min_edge<<", label = "<<treeT.treeNodes[i].left<<", fontsize =30]";					
				
				if(treeT.isTip(treeT.treeNodes[i].left))
				{      ss<<" n_"<<i<<" -- "<<treeT.lineageNodeHash[treeT.treeNodes[i].left]+1<<"}"<<endl;}
				else
				{      ss<<" n_"<<i<<" -- n_"<<treeT.treeNodes[i].left<<"}"<<endl;}
			}		
		}
	}
	ss<<endl<<"}"<<endl;
	
	return ss.str();
}

double lineage::logBinomial(long i, long j)
{
	return gsl_sf_lnchoose((unsigned int)i, (unsigned int)j);
}

double lineage::logFactorial(long n)
{
	if ((n==1)||(n==0))
	{ 	return 0;	}
	else
	{  	return gsl_sf_lnfact((unsigned int)n);	}
}

double lineage::glog(double value)
{
	if (value > 1e-323)
	{
		return(gsl_sf_log(value));
	}
	else
	{
		return(-743.75);
	}
}
#endif /* LINEAGE_H */
